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A  PARAMETRIC  PATCH  SURFACE  GEOMETRY 
DEFINITION  FOR  THE  THREE-DIMENSIONAL 
RETARDED  POTENTIAL  TECHNIQUE 


INTRODUCTION 

This  study  is  a  part  of  the  effort  to  apply  the  retarded  potential  integral  to  three-dimensional  prob¬ 
lems  in  fluid-structure  interaction  for  arbitrary  bounding  geometries.  This  paper  will  be  confined  to  a 
discussion  of  the  geometric  aspects  of  retarded  potential  computations.  A  related  paper  covers  the 
application  of  the  three  dimensional  technique  to  problems  of  shock  [61. 

Mitzner  [1.2]  developed  a  method  for  discretizing  the  integral  spatially  and  temporally,  and 
applied  the  technique  to  the  problem  of  a  spherical  rigid  submerged  body  with  two-dimensional  axisym- 
metry.  Huang.  Everstine.  and  Wang  [3]  extended  the  method  to  handle  a  non-rigid  boundary  condi¬ 
tion.  and  tested  the  case  of  a  spherical  elastic  shell,  modeled  axisymmetrically. 

This  work  builds  upon  the  computer  program  developed  by  Mitzner  and  extended  to  non-rigid 
boundaries  by  Huang,  Everstine,  and  Wang.  A  fully  three-  dimensional  boundary  element  formulation 
for  arbitrarily  shaped  surfaces  has  been  implemented.  The  entire  approach  is  designed  to  be  compatible 
with  finite  element  methods,  important  in  the  case  of  coupled  problems  involving  non-rigid  submerged 
structures. 

In  the  new  development,  the  surface  geometry  is  totally  defined  by  the  boundary  element  input 
data.  The  Cartesian  coordinates,  curvatures,  and  surface  coordinate  base  vectors,  which  are  needed  at 
various  points  by  the  coefficients  calculation  routine,  are  all  derived  from  the  natural  surface  coordinate 
interpolation  functions  of  the  boundary  patch.  The  numerical  definition  of  geometry  is  required  for  a 
problem- independent  retarded  potential  analysis  capability,  particularly  when  coupled  to  a  finite  element 
program.  This  avoids  the  need  to  individually  tailor  surface  parametric  equations  for  each  problem 
geometry  encountered,  which  would  require  reprogramming.  Beyond  this  consideration,  for  an  arbi¬ 
trary  geometry  it  may  be  difficult  to  obtain  a  set  of  closed-form  parametric  equations  to  define  a  partic¬ 
ular  problem 

RETARDED  POTENTIAL  FORMULATION 

A  discussion  of  the  retarded  potential  equation  and  its  discretized  approximation  may  be  found  in 
Refs.  [1.2.31  with  a  derivation  in  Ref.[4l.  The  form  of  the  equation  for  calculation  of  total  pressures, 
p.  on  the  boundary  of  a  smooth,  closed  surface  in  an  infinite  fluid  is  as  follows: 

plx.i)  -  2p"r<I.r)  -  p/2w  /  wix  -n  js- 
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where  p  and  c  are  the  fluid  density  and  sonic  velocity,  respectively,  x  is  the  surface  position  of 
the  field  point,  x'  is  the  surface  position  of  an  integration  point,  R  -  lx  —  x'l  is  the  distance  between 

D 

field  and  integration  points,  t‘  -  t - is  the  retarted  time  of  the  integration  point,  pm  is  the  incident 

c 

pressure  field  in  the  fluid,  n‘  is  the  surface  normal  at  the  integration  point  directed  out  of  the  fluid,  and 
ii'  is  the  normal  direction  acceleration  of  the  surface. 


From  the  integral,  the  basic  geometric  influence  factors  to  be  evaluated  are  as  follows: 


R  - 

and 


•  1  x  —  xl 
dS'/R 

R2  9/i 


(3a) 

(3b) 

(3c) 


Computation  of  the  discretized  forms  of  these  for  the  parametric  surface  patch  is  the  subject  of  this 
report. 


DISCRETIZATION 

For  computation,  the  retarded  potential  integral  is  discretized  into  a  set  of  algebraic  equations. 
The  surface  pressure  field  is  approximated  by  subdividing  the  surface  into  zones  of  spatially  constant 
pressure.  The  surface  normal  acceleration  field,  if  not  identically  zero,  is  also  approximated  as  constant 
on  the  surface  of  the  same  zone.  Once  the  pressure  and  acceleration  (and  pressure  derivative)  are 
approximated  as  constant  over  a  part  of  the  surface,  these  quantities  can  be  factored  out  of  the  integral 
expressions.  This  results  in  integrals  which  are  dependent  on  surface  geometry  alone.  There  is  a  need 
to  further  subdivide  the  surface  of  zones  into  a  mesh  of  subzone  elements  in  order  to  accurately 
numerically  integrate  the  geometric  influence  factors,  as  well  as  to  accurately  obtain  the  time-delayed 
influences  of  pressures  and  accelerations  from  each  preceding  time  step  [1.2.31.  Then,  the  acceleration 
influence  and  pressure  influence  integrals  over  the  subzone  element.  5,.  are  as  follows: 

{  f 

f  dS*  <4b> 

i  R2  9h 

Discussion  of  the  discretization  of  retarded  time  into  steps  is  left  to  the  related  report  on  the  application 
of  the  technique  to  shock  response  problems  (61  or  can  be  found  in  the  original  references  [1.2.31 


In  the  original  solution  procedure,  the  computations  are  based  on  a  surface  coordinate  system. 
(a,jS),  common  to  the  entire  boundary  and,  unique  to  the  to  the  pirticular  geometry  Then  the  subdi¬ 
vision  of  the  surface  into  zones  and  subzone  elements  is  made  to  follow  the  lines  of  the  coordinate  sys¬ 
tem  (1,2).  The  flexibility  of  the  new  development  allows  for  a  series  of  independent  surface  coordinate 
systems,  one  per  surface  patch,  as  will  be  described  later  In  the  following  discussion  the  dimensions  of 
a  subzone  element,  expressed  in  surface  coordinates,  are  la  by  10.  and  the  element's  surface  area  is 

V 


Computationally,  the  subzone  elements  are  classified  into  two  types,  singular  and  non-singular 
The  area  of  a  non-singular  element  is  approximated  as  a  flat  surface  ungem  to  the  integration  point  ji 
its  center.  With  the  boundaries  of  the  element  following  the  lines  of  the  surface  coordinate  system,  the 
cross  product  of  the  Euclidean  base  vectors  becomes  the  key  quantity  in  obuming  the  physical  size  of 
the  element  and  iu  normal-direction  vector  The  resulting  expressions  for  J  JS  R  and 

/  1/*J  9R/9n  dS 
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'  over  the  non-singular  element  are 


and 
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respectively.  Distance  measurement  between  a  non-singular  element  and  the  field  point  of  a  zone  is 
approximated  by  locating  the  field  point  at  the  center  of  the  field  zone  and  the  integration  point  at  the 
center  of  the  integration  subzone  element. 

For  the  surface  area  immediately  surrounding  the  field  point,  modified  calculations  must  be  used 
to  compute  the  pressure  and  acceleration  integration  influences,  because  of  the  close  proximity  of  the 
field  point  and  the  presence  of  a  singularity  at  the  field  point  itself.  For  the  pressure  integral,  the  sur¬ 
face  surrounding  the  field  point  is  approximated  by  an  osculating  paraboloid,  with  principal  curvatures 
matching  at  the  central  field  point,  and  an  exact  integral  is  evaluated  over  this  surface.  In  the  current 
implementation,  this  integral  has  been  simplified  by  requiring  the  surface  coordinate  system  to  be 
aligned  with  the  principal  curvature  directions.  With  the  simplification,  the  expression  for 
f  l/R2  dR/dn'  dS'  over  the  singular  element  becomes 
5, 

2UaYjn^^-  +  K9X0ln-~r-2-)  111  (6) 

* 0  * 0 

where 

Xt  -  1/2  IjJ  Aa 
Y ,  -  1/2  l?J  A/3 

and 

-  (Xf+Y2)* 

In  the  above  formulas  the  <c’s  are  the  principal  curvatures  at  the  field  zone  center,  with  the  rest  of  the 
quantities  derived  from  the  Euclidean  base  vectors  and  the  element  size,  as  measured  in  surface  coordi¬ 
nates.  The  singular  element  computations  require  the  distance,  R,  in  the  expression  R/c  of  the  pres¬ 
sure  derivative  term.  In  this  case  the  integration  point  is  taken  to  be  a  distance  of  one-fourth  the  diag¬ 
onal  dimension  of  the  singular  subzone  from  the  zone-centered  field  point  i'/ir0). 

For  the  acceleration  influence,  a  cruder  approximate  integral  is  used  here.  This  is  acceptable 
because  the  weaker  singularity  yields  a  less  significant  influence  of  the  near-by  accelerations.  The 
approximation  involves  replacing  the  rectangular  element  with  a  circle  of  the  same  area,  and  integrating 
the  expression  dS’/ R  over  a  flat  surface.  The  result  is  MwX,  Y,)'H  [3al.  As  can  be  seen  from  the  above 
expressions,  the  influence  coefficients  calculations  require  the  evaluation  of  three  basic  local  geometric 
quantities:  the  position  vector,  the  coordinate  base  vectors,  and  the  curvatures.  The  subzone  element 
sues  are  determined  from  the  user— specified  subzone  element  mesh  refinement. 

PARAMETRIC  PATCH  SURFACE  DEFINITION 

The  surface  geometry  definition  implemented  for  the  retarded  potential  method  is  the  same  as 
that  used  for  isoparametric  finite  element  shells  in  three  dimensional  space  (5)  The  boundary  surface 
is  divided  into  patches  to  coincide  with  zones  of  constant  pressure.  When  coupled  to  a  structural 
analysis  code  .these  also  coincide  with  the  shell  surface  finite  elements,  allowing  one  surface  definition 
to  be  given  to  both  the  structural  problem  and  the  fluid  pressure  problem  Each  surface  patch  has  its 
own  curvilinear,  non-dimensional  coordinate  system,  which  follows  the  surface  contour  with  coordinate 
values  varying  from  -1  to  -H  This  is  referred  to  as  a  natural  coordinate  system,  with  coordinate  direc¬ 
tions  r  and  s.  The  surface  zone  used  in  this  study  has  a  3  *  3  arrangement  of  nodal  points  to  match 
the  9-noded  version  of  the  ADINA  shell  finite  element.  The  nodes  lie  at  whole  number  values  of  (he 
parametric  coordinate  pairs,  designated  <  /■ . sf )  for  node  /  The  nodal  coordinates.  ?v.  define  the  Carte¬ 
sian  coordinates  of  the  entire  element  surface  through  biquadratic  interpolation  functions  of  the  two 
natural  coordinates,  as  shown  in  Eq  (7) 
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(7) 


X(r,s)  -  £//,(r.s)  Xf 

i 

The  interpolation  functions  are  built  up  from  elementary  terms  of  the  form,  F,  shown  as  Eq.  (8), 
which,  in  turn,  are  products  of  the  unidirectional  factors,  Q,  shown  in  Eqs.  (8a). 

F(r,s,ri,Si )  -  Q(/’,r/)  Q(s,s,)  (8) 

for  |r,|-»l  Q(.rtr ,)  —  1 h 

for  Ir,!— 0  Q(r,rt)  -  (1— r2)  (8a) 

The  function,  F,  has  a  value  of  one  at  the  node  i  and  a  value  of  zero  at  all  edges  of  the  finite  element 
that  do  not  touch  node  i.  For  edge  and  corner  nodes,  achieving  an  interpolation  function  that  has  a 
zero  value  at  all  nodes  other  than  /  requires  subtraction  of  fractional  values  of  other  interpolation  func¬ 
tions  from  the  F  expression.  The  full  interpolation  functions  are  related  to  the  F  expressions  and  to 
each  other  by  the  following  expressions. 

Hi  -  F(.r,s,  1,  1)  -  .5  (Hs+Ht)  -  .25  H9 

H2  -  F(r,s,- 1,  1)  -  .5  (Hs+H6)  -  .25  H9 

Hi  -  F(r,s,- 1,-1)  -  .5  (tf4+//7)  -  .25  H9 

Ha  -  F(r,s,  1,-1)  -  .5  (//7+//8)  -  .25  H9 

Hs  -  F(r,s,  0,  1)  -  .5  H9  (9) 

H6-  F(r,j,-1,0)  -.5  H9 

Ht- F(r,s,  0,-1)  -.5  H9 

Ht  -  F(r,s ,  1,  0)  -  .5  H9 

H9-  F(r,s,  1,-1) 
where  H,  -  H(r,s). 


GEOMETRIC  COMPUTATIONS 


As  discussed  earlier,  three  types  of  geometric  information  are  needed  for  the  calculation  of  the 
retarded  potential  influence  coefficients. 

1.  The  Cartesian  Position  Vectors  are  needed  to  locate  the  field  and  integration  points  in  space. 

2.  Base  Vectors  (which  the  derivatives  of  the  position  vector  with  respect  to  the  natural  coordinate 
directions)  are  needed  to  determine  the  orientation  and  size  of  the  integration  subzones. 


3.  Principal  Surface  Curvatures  are  required  at  the  field  points  to  determine  the  osculating  paraboloid 
for  the  singular  subzone  calculation  discussed  earlier. 

The  expressions  for  the  base  vectors,  g,  and  g, ,  can  be  obtained  by  taking  the  r -  and  s-direction 
derivatives  (i.e.  a  derivative  with  respect  to  one  of  the  natural  coordinates)  of  the  position  vector  inter¬ 
polation  function  and  are  given  by: 


where  q  is  r  or  j. 


The  rate  of  change  in  surface  orientation  is  expressed  as  the  rate  of  change  of  a  unit  normal  vec¬ 
tor  with  respect  to  a  position  change  in  a  given  direction  along  the  surface.  Because  the  unit  vector  can 
only  change  in  direction  and  not  length,  its  change  from  one  point  to  another  is  a  measure  of  angle 
change,  with  the  added  feature  of  showing  the  direction  of  the  angle  change.  The  component  of  the 

angle  change  in  the  direction  of  position  change  gives  the  angle  of  deflection.  If  one  considers  a  rate  of 

change  in  the  unit  normal  instead,  then  the  component  in  the  direction  of  position  change  becomes  the 
curvature.  This  component  is  obtained  (in  Eq.  (ID)  by  taking  the  dot  product  of  the  rate  of  change  of 
the  unit  vector  and  the  unit  vector  in  the  direction  of  travel.  [7] 

dq  -  nn 

*<>  ~  — T~  '  '»  (D) 

"  ds„  ' 

where  q  -  either  one  of  the  natural  coordinate  directions,  r  or  s 
Sq  —  physical  distance  in  the  q  direction 
i„  —  unit  vector  normal  to  the  surface 
/,  —  unit  vector  in  the  q  direction 

The  derivation  of  the  factors  in  Eq.  (11)  are  as  follows: 

—  — i-7  where  g.  -  the  base  vector  in  the  q  direction 

If,  I 


/,  gq/  '  gq  I 


Z  -  vj\v. 


where 


f,  *  f. 


"wn  /i—  |  _  ,  -  ■* 

_  .  _/UJ  . 


where 


Bgr  - 

-a7x*  • 


af,  - 
g' 
Bq 


As  can  be  seen,  the  curvature  calculation  requires  the  evaluation  of  the  the  base  vector  derivatives,  i.e. 
the  second  derivatives  of  the  position  vector  with  respect  to  the  surface  coordinates. 

All  the  calculations  for  position  vector,  base  vectors,  and  curvatures  can  use  the  same  interpola¬ 
tion  scheme  expressed  by  Equations  (S).  By  simply  taking  the  appropriate  derivatives  of  the  F  terms, 
the  corresponding  appropriate  derivatives  of  the  H  functions  are  obtained.  These,  in  turn,  make  up  the 
expressions  for  all  the  quantities  needed. 

TEST  PROBLEM 

The  retarded  potential  technique  was  studied  for  the  rigid  scattering  problem  to  verify  the  3- 
dimensional  integration  method.  The  test  was  run  to  duplicate  the  Gaussian  pulse  loading  on  a  sphere 
case  that  Mitzner  [1]  tested.  The  results  obtained,  plotted  in  Fig.  3,  are  virtually  identical  to  the  separa¬ 
tion  of  variables  solution,  validating  the  numerical  geometry  formulation. 


The  structural  model  used  for  the  rigid  body  test  problem  is  shown  in  Fig.  1.  This  is  a  quarter 
section  of  a  sphere  used  to  model  an  entire  spherical  rigid  body  with  symmetry  conditions  imposed  on 
the  two  cutting  planes.  The  retarded  potential  program  allows  for  such  symmetry  conditions.  The  sym¬ 
metries  are  imposed  to  save  computation  cost  and  can  be  allowed  because  both  the  test  structure  and 
the  test  loadings  are  axisymmetric  with  respect  to  the  loading  direction.  The  subdivision  into  pressure 


zones  is  more  refined  in  the  direction  of  wave  travel  in  order  to  capture  the  detailed  response  expected 
in  that  direction.  Subzone  division  for  this  problem  is  two  refinements,  one  for  the  influence  between 
two  different  pressure  zones  and  one  for  the  influence  of  a  pressure  zone  on  itself.  These  are  shown  in 
Figs.  2. 


The  computed  curvat  res  at  the  centers  of  the  pressure  zones 
are  as  follows: 
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9.0 

0.01027 

0.010043 

22.5 

0.01040 

0.010015 

31.5 

0.01040 

0.010015 

40.5 

0.01040 

0.010015 

49.5 

0.01040 

0.010015 

58.5 

0.01040 

0.010016 

67.5 

0.01040 

0.010015 

76.5 

0.01040 

0.010016 

85.5 

0.01040 

0.010015 

The  list  of  k4  shows  that,  even  with  parametric  surfaces  stretched  over  45  degrees  of  circumference, 
the  curvature  values  are  reasonably  close  to  the  exact  value,  0.01,  for  a  sphere  of  radius  “  100.,  and 
the  values  of  k«  are  virtually  exact. 

CONCLUSIONS 

This  study  has  developed  a  parametric  patch  formulation  for  the  surface  geometric  computations 
of  the  discretized  retarded  potential  integral.  The  parametric  patch  definition  establishes  the  retarded 
potential  integral  as  a  fully  developed  boundary  element  formulation.  With  the  added  programming 
enhancement  to  eliminate  the  requirement  that  surface  patch  grid  lines  follow  the  principal  curvature 
directions  on  the  surface,  complete  flexibility  in  surface  geometry  specification  will  be  achieved. 
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